%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Light-shift
% Stationary optical Bloch equations for a Lambda system,
%  n = [22 33 12 21 13 31 23 32],
%	   --- |2> 
%	   / \			
%	  /  --- |3>
%	 /	  		  
%  --- |1> 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LightShift
    figure(1);          clf;
    col='kbgrmcy';

    Gam=2;
    omg=[0 .49 .6 .8 1]*Gam;
    Del=[-1:.01:1]'*Gam;

    % numerical calculation

    for kk=1:length(omg)
        Omg=omg(kk);
        for nn=1:length(Del)
            Eigen = eig( [0 Omg/2 ; Omg/2 Del(nn)-1i*Gam/2] );
            [Shift(:,nn),I] = sort(real(Eigen));
            Width(:,nn) = -imag(Eigen(I));
            if kk<3
                Shift(:,nn) = real(Eigen);
                Width(:,nn) = -imag(Eigen);
            end
        end
        E1(:,kk)=Shift(1,:);        W1(:,kk)=Width(1,:);
        E2(:,kk)=Shift(2,:);        W2(:,kk)=Width(2,:);
    end

    [omg2,Del2]=meshgrid(omg,Del);
    subplot(221),
    for kk=1:length(omg)
        plot(Del,E1(:,kk),col(kk),Del,E2(:,kk),col(kk),'LineWidth',1)
        hold on;
    end
    hold on;
    set(gca,'xLim',[-1 1]*Gam,'yLim',[-1 1]*Gam,'FontSize',16,'LineWidth',1)
    xlabel('\fontsize{16}\Delta/\Gamma')
    ylabel('\fontsize{16}\Delta{\itE}/\Gamma')
    text(-3,2,'\fontsize{16}(a)')
    subplot(222),
    for kk=1:length(omg)
        plot(Del,W1(:,kk),col(kk),Del,W2(:,kk),col(kk),'LineWidth',1)
        hold on;
    end
    set(gca,'FontSize',16,'LineWidth',1)
    xlabel('\fontsize{16}\Delta/\Gamma')
    ylabel('\fontsize{16}\Delta\Gamma')
    text(-3,1,'\fontsize{16}(b)')

    % exact analytical calculation (are identical to numerical)

    Lmb=Del2-1i*Gam/2;
    GG =sqrt(omg2.^2+Lmb.^2);
    Shift1p=Lmb/2+GG/2;
    Shift1m=Lmb/2-GG/2;
    subplot(221),
    for kk=1:length(omg)
    %      plot(Del2,real(Shift1p(:,kk)),col(kk),Del2,real(Shift1m(:,kk)),col(kk),'LineWidth',1)
    end
    G  =sqrt(omg2.^2+Del2.^2);
    Shift2p=1/2*Del2+1/2*sqrt(sqrt(G.^4/4+Gam^2/8.*(Del2.^2+Gam^2/8-omg2.^2))+G.^2/2-Gam^2/8);
    Shift2m=1/2*Del2-1/2*sqrt(sqrt(G.^4/4+Gam^2/8.*(G.^2+Gam^2/8-2*omg2.^2))+G.^2/2-Gam^2/8);

    % weak coupling aproximation

    s  =2*omg2.^2./(4*Del2.^2+Gam^2);
    Shift3p=Lmb/2+Lmb/2.*(1+s);                                                
    Shift3m=Lmb/2-Lmb/2.*(1+s);

    % strong coupling aproximation

     del=1/2*Del2+sqrt(Del2.^2+Gam^2);                                          
    G1 =1i/2*Gam*(1-Del2./GG);
    G2 =1i/2*Gam*(1+Del2./GG);
    Shift4p=1/2*Del2+1/2*G-Gam^2*omg2.^2/16./G.^3;
    Shift4m=1/2*Del2-1/2*G+Gam^2*omg2.^2/16./G.^3;

    subplot(221),
    for kk=1:length(omg)
    %     plot(Del2,real(Shift4p(:,kk)),['.',col(kk)],Del2,real(Shift4m(:,kk)),['.',col(kk)],'LineWidth',1)
    end

    Bloch(.4,3);
    Bloch(2,4);
end


function Bloch(o12,pic)
    G12=              1;                                                    % decay of populations
    G23=              1;
    G13=	        .01;
%     o12=          6*G12;													% Rabi frequencies
    o23=	      .2*G23;
    D12=[-1:.01:1]*G12*2;                                                   % detunings
    D23=[-1:.01:1]*G12*2;
    d23=	           0;

    for kk=1:length(D23)
      d23=	 D23(kk);
      for ll=1:length(D12)
            d12= D12(ll);
            n  = Matrix_L(G12,G23,G13,o12,o23,d12,d23);
            nt(kk,ll)=real(n(1));
        end
    end

    subplot(2,2,pic),surf(D23,D12,nt');
    set(gca,'FontSize',16,'LineWidth',1.3)
    xlabel('\fontsize{16}\Delta_{23}/\Gamma_{23}');
    ylabel('\fontsize{16}\Delta_{12}/\Gamma_{12}');
    zlabel('\fontsize{16}{\it\rho_{22}}');
    if pic==3
        txt='(c)';
    else
        txt='(d)';
    end        
    text(2.5,-3,0,['\fontsize{16}',txt])
    shading flat
    view([90 -90])
    axis tight;
end


function n = Matrix_L(G12,G23,G13,o12,o23,d12,d23)
    I=          1i/2;
    g12=     G12+G23;                                                       % decay of coherences
    g13=         G13;
    g23= G12+G23+G13;
    d13=     d12-d23;
    L12=g12/2+1i*d12;
    L13=g13/2+1i*d13;
    L23=g23/2+1i*d23;

    b=[0 0 -I*o12 I*o12 0 0 0 0]';
    
    M=[
       -G12-G23	0		-I*o12	I*o12	0		0		I*o23	-I*o23
        G23		-G13	0		0		0		0		-I*o23	I*o23
        -1i*o12	-I*o12	-L12	0		I*o23	0		0		0
        1i*o12	I*o12	0		-L12'	0		-I*o23	0		0
        0		0		I*o23	0		-L13	0		-I*o12	0
        0		0		0		-I*o23	0		-L13'	0		I*o12
        I*o23	-I*o23	0		0		-I*o12	0		-L23'	0
        -I*o23	I*o23	0		0		0		I*o12	0		-L23
    ];

    n = -M\b;
end